library(Rcpp)
library(RcppEigen)
install.packages("npideal", repos="https://r.tahk.us/")
## To install from source:
## devtools::install_github("atahk/npideal", ref="v0.1.1")
library(npideal)

library(pscl)
library(xtable)
library(parallel)

num.vote.groups <- c(100,500,1000)
simfile <- "mc-2-out.RData"

sim.data <- function(prob.fun, k, n, pxf) {
    p1 <- t(sapply(rep(k,n),pxf))
    ## polarity <- sapply(p1[,ncol(p1)]-p1[,1],sign)
    ## mean.prob <- apply(p1*polarity, 2, mean)
    v1 <- array(rbinom(length(p1),1,p1),dim=dim(p1))*2-1
    ## attr(v1, "mean.prob") <- mean.prob
    return(v1)
}

pxf.sym <- function(x, sample.polarity=TRUE)
    pnorm(1:x,sample(x,1),2,lower.tail=sample(c(!sample.polarity,TRUE),1))
      
pxf.asym <- function(x, sample.polarity=TRUE)
    punif(.5^(1:x), lower.tail=sample(c(!sample.polarity,TRUE),1))

sim.model <- function(num.votes, k, pxf, swap.fun) {
    num.votes1 <- num.votes2 <- num.votes
    v1 <- sim.data(pxf, k, num.votes1, pxf)
    v2 <- sim.data(pxf, k, num.votes2, pxf)
    if (!is.null(swap.fun)) {
        swap.which <- swap.fun(k)
        v2[,sort(swap.which)] <- v2[,swap.which]
    }
    return(votetest(v1,v2,as.integer(NULL),numLegisShuffles=3000,numBootstrapSims=200))
}


swap.fun2 <- function(x) return(sample(x-1,1)+(1:0))
swap.fun3 <- function(x) sample(x)

getOption("mc.cores")
num.sim <- 1000
new.est <- FALSE
if (file.exists(simfile)) {
    load(simfile)
} else {
    if (file.exists(paste0(simfile,".tmp"))) {
        load(paste0(simfile,".tmp"))
    } else {
        models <- list()
    }
    for (num.votes in num.vote.groups) {
        modname <- paste("None",num.votes,sep="-")
        if (!(modname %in% names(models))) {
            models[[modname]] <- mclapply(rep(num.votes,num.sim), sim.model, k=9, pxf=pxf.sym, swap.fun=NULL)
            new.est <- TRUE
            save(models, file=paste0(simfile,".tmp"), compress=TRUE)
        }
        modname <- paste("Min",num.votes,sep="-")
        if (!(modname %in% names(models))) {
            models[[modname]] <- mclapply(rep(num.votes,num.sim), sim.model, k=9, pxf=pxf.sym, swap.fun=swap.fun2)
            new.est <- TRUE
            save(models, file=paste0(simfile,".tmp"), compress=TRUE)
        }
        modname <- paste("Max",num.votes,sep="-")
        if (!(modname %in% names(models))) {
            models[[modname]] <- mclapply(rep(num.votes,num.sim), sim.model, k=9, pxf=pxf.sym, swap.fun=swap.fun3)
            new.est <- TRUE
            save(models, file=paste0(simfile,".tmp"), compress=TRUE)
        }
    }
    if (new.est) {
        save(models, file=simfile, compress=TRUE)
        unlink(paste0(simfile,".tmp"))
    }
}

texexpr <- function(x)
    as.expression(as.formula(paste0(gsub(" ", "~", gsub("\\\\emph\\{(.*?)\\}", "italic(\\1)", x)),"~a"))[[2]])

pdf("pvals-qqplot.pdf", width=10, height=10)

layout(cbind(c(0,1:3),rbind(4:6,matrix(7:15,3,3))), widths=c(2,3,3,3), heights=c(1,3,3,3))
par(mar=c(0,0,0,0))
cex.lab <- 1.5
for (num.votes in num.vote.groups) {
    plot.new()
    text(.5,.5,paste(num.votes,"votes\nper group"),cex=cex.lab)
}
plot.new()
text(.5,.5,paste("No difference\nbetween groups"),cex=cex.lab)
plot.new()
text(.5,.5,paste("Minimal differences\n(two neighboring\nlegislators switch order)"),cex=cex.lab)
plot.new()
text(.5,.5,paste("Large differences\n(independent orders\nfor each group)"),cex=cex.lab)
## par(mar=c(4,2.5,4,2.5))
par(mar=c(4,4,1,1))
for (change in c("None","Min","Max")) {
    for (num.votes in num.vote.groups) {
        modname <- paste(change,num.votes,sep="-")
        model <- models[[modname]]
        pval <- sapply(model, function(x) x$p)
        ## pval <- sapply(model, function(x) pchisq(x$dev,2*x$df,lower.tail=FALSE))
        qqplot(qunif(ppoints(length(pval))), pval,
               xlab = "Theoretical quantiles of uniform dist.",
               ylab = texexpr("Empirical quantiles of \\emph{p}-values"),
               xlim = c(0,1), ylim = c(0,1), cex=.75, pch=16)
        abline(c(0,1), col="red")        
        ## hist(pval, main='', breaks=seq(0,1,.05),
        ##      col="grey",yaxt='n',xaxs='i',yaxs='i',
        ##      xlab="\\emph{p}-value")## xlab=expression(italic(p)-value))
        ## ## plot(density(pval,from=0,to=1),bty='n',yaxt='n',ylab='')
    }
}

dev.off()


pdf("pvals.pdf", width=8, height=8)

layout(cbind(c(0,1:3),rbind(4:6,matrix(7:15,3,3))), widths=c(2,3,3,3), heights=c(1,3,3,3))
par(mar=c(0,0,0,0))
cex.lab <- 1.5
for (num.votes in num.vote.groups) {
    plot.new()
    text(.5,.5,paste(num.votes,"votes\nper group"),cex=cex.lab)
}
plot.new()
text(.5,.5,paste("No difference\nbetween groups"),cex=cex.lab)
plot.new()
text(.5,.5,paste("Minimal differences\n(two neighboring\nlegislators switch order)"),cex=cex.lab)
plot.new()
text(.5,.5,paste("Large differences\n(independent orders\nfor each group)"),cex=cex.lab)
## par(mar=c(4,2.5,4,2.5))
par(mar=c(4,4,1,1))
for (change in c("None","Min","Max")) {
    for (num.votes in num.vote.groups) {
        modname <- paste(change,num.votes,sep="-")
        model <- models[[modname]]
        pval <- sapply(model, function(x) x$p)
        hist(pval, main='', breaks=seq(0,1,.05),
             col="grey",yaxt='n',xaxs='i',yaxs='i',ylab="",
             xlab=texexpr("\\emph{p}-value"))
        ## plot(density(pval,from=0,to=1),bty='n',yaxt='n',ylab='')
    }
}

dev.off()


models.p <- lapply(models, function(x) sapply(x, function(y) y$p))

alpha <- .05
num.sims <- length(models.p[[1]])
simvotes <- 
list(simvotesone = num.vote.groups[1],
     simvotestwo = num.vote.groups[2],
     simvotesthree = num.vote.groups[3])
simprobs <-
list(simvotesonecoltwoprob =
         mean(models.p[[paste("Min",num.vote.groups[1],sep="-")]]<alpha),
     simvotestwocoltwoprob =
         mean(models.p[[paste("Min",num.vote.groups[2],sep="-")]]<alpha),
     simvotesthreecoltwoprob =
         mean(models.p[[paste("Min",num.vote.groups[3],sep="-")]]<alpha),
     simvotesonecolthreeprob =
         mean(models.p[[paste("Max",num.vote.groups[1],sep="-")]]<alpha),
     simvotestwocolthreeprob =
         mean(models.p[[paste("Max",num.vote.groups[2],sep="-")]]<alpha),
     simvotesthreecolthreeprob =
         mean(models.p[[paste("Max",num.vote.groups[3],sep="-")]]<alpha))

siminfo <- "siminfo.tex"
cat("\\newcommand\\alphapct{", 100*alpha, "}\n", sep="", file=siminfo)
for (i in names(simvotes)) {
    cat("\\newcommand\\",i,"{",simvotes[[i]],"}\n",sep="",file=siminfo,append=TRUE)
}
for (i in names(simprobs)) {
    probs <- simprobs[[i]]
    if (probs == 1)
        fout <-
            sprintf("over %.1f", 100*(num.sims-1)/num.sims)
    else
        fout <-
            sprintf("%.1f", 100*probs)
    cat("\\newcommand\\",i,"{",fout,"}\n",sep="",file=siminfo,append=TRUE)
}
